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Abstract 

Mathematical physics problems are often formulated using differential 
oprators of vector analysis - invariant operators of first order, namely, di- 
vergence, gradient and rotor operators. In approximate solution of such 
problems it is natural to employ similar operator formulations for grid prob- 
lems, too. The VAGO (Vector Analysis Grid Operators) method is based on 
such a methodology. In this paper the vector analysis difference operators 
are constructed using the Delaunay triangulation and the Voronoi diagrams. 
Further the VAGO method is used to solve approximately boundary value 
problems for the general elliptic equation of second order. In the convection- 
diffusion-reaction equation the diffusion coefficient is a symmetric tensor of 
second order. 

Keywords: finite difference method, unstructured grids, Delaunay 
triangulation, Voronoi diagrams, convection-diffusion problems. 
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1. Introduction 

In the theory of difference schemes [l| the problem of approximation of 
a differential problem via grid one received primary emphasis. A universal 
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method of balance (integro-interpolation method) is used as the basic one to 
construct discrete analogs. It was proposed bySamarskii ^ and currently 
it is known as the control volume method js], 0, Isf. A difference scheme is 
constructed using the balance method via integration of the initial equation 
over a control volume — a part of the computational domain adjacent to the 
specified grid point. 

In solving mathematical physics problems on general unstructured grids 
there is often used the Delaunay triangulation which displays some optimal 
features [6|. For the Delaunay triangulation it is natural to consider as control 
volumes a Voronoi polygon — a set of points lying closer to the considered 
grid point in compare with all others. Consistent approximations of the 
vector analysis operators on the basis of the Delaunay triangulation were 
constructed in [3] for some problems of mathematical physics. A systematic 
investigation of approxinations for gradient and divergence operators on dual 
Voronoi-Delaunay partitionings has been conducted in j^, 0] for model 2D 
and 3D problems. 

The present work continues our studies on Vector Analysis Grid Oper- 



ators method. In paper [lOj there are discussed problems of constructing 



consistent approximations for basic operators of vector analysis — gradient, 
divergence amd rotor operators. Special representations were derived for 
truncation errors that allowed to obtain the corresponding estimates for the 
convergence rate of the approximate solution to exact one. Possibilities of the 
VAGO method have been demonstrated on some scalar and vector problems. 
Difference schemes for steady-state convection-diffusion problems were con- 



structed on the basis of the dual Voronoi-Delaunay partitioning in [ll|, |12 



In these problems (see, e.g., [l3j) the emphasis is on properties of discrete op- 



erators for the convection and diffusion transport. The convective terms are 
written lj,ll2| in the divergent (conservative), non-divergent (characterictic) 



and symmetric forms. In \XQl there were designed VAGO approximations as 
well as investigated truncation errors for the Dirichlet problem for stationary 
and usteady convection-diffusion problems for isotropic media. 

Paricular emphasis should be placed on approximation of diffusion trans- 
port operator for an anisotropic medium where the diffusion coefficient is 
a symmetric tensor of second order. In this case we should control at the 
discrete level not only the seld-adjoint property of the grid diffusion opera- 
tor but the monotonicity and conservation properties for the corresponding 
grid approximatioms, too [T]. The present-day view on approximation of 
problems with mixed derivatives on rectangular grids is given in [l5|. Prop- 
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erties of standard finite-element approximations on triangular or tetrahedral 
meshes are discussed in 16|] in detail. To obtain the monotonicity property, a 
modification of standard linear Galerkin finite element discretizations of the 
Laplace operator was performed. Much attention is given to constructing 
control volume schemes for the approximate solution of convection-diffusion 



boundary value problems with an anisotropic diffusion coefficient (see [17 
and references therein). In particular (see also [ISj]), the control volume 
schemes are developed with an anisotropic diffusion for approximating a lo- 
cal discrete gradient. The results of work [8j] are generalized in [1^ to the 
2D div-curl system in anisotropic media. 

The outline of this paper is the following. In Section 2, the boundary- 
value problem for the general elliptic equation of second order is formulated 
in an invariant form as the problem for a convection-diffusion equation based 
on using vector analysis operators — divergence and gradient. The diffusion 
coefficient is a symmetric positive definite tensor of second order. On dual 
Voronoi-Delaunay partitionings (Section 3) there are introduced the corre- 
sponding spaces of grid functions. Using the VAGO method the difference 
operators of divergence and gradient are constracted and investigated in Sec- 
tion 4. Difference schemes for the considered convection-diffusion-reaction 
problem are studied in Section 5. Approximation of the diffusion transport 
operator with an anisotropic diffusion coefficient has reseived much consid- 
eration. 



2. The boundary value problem 

The boundary value problem is considered for an elliptic equation of sec- 
ond order. In n-dimensional bounded domain Q with smooth enough bound- 
ary dQ function m(x), x = (xi, x„) satisfies the following equation 



a, 3=1 a=l a=l 



du ^ 



c(x)n(x) = /(x), xGl]. (1) 
The coefficients at the higher order derivatives meet the ellipticity condition: 



a, 13=1 a=l 



n 

2 
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6>0, VyGR", VxG^]. (2) 
For other coefficients of equation ([T]) we have 

iE%^ + cW>0. (3) 

Assume that the right hand side of the equation can be represented as follows 

/(x) = /o(x) + X:f^. (4) 

The simplest boundary value problem is considered for ([1]) with the following 
boundary conditions 

u(x) = 0, X G d^. (5) 

In Hilbert space L^i^C) the scalar product and norm are introduced like 
this 

[u,v) = \ u(x)t;(x)dx, ll-ull = (n, n)"*"/^. 



In 7/(grad) we introduce 

(m, ^^)w(grad) = {n^V) + ^ , ||w||w(grad) = ^(grad)) • 

For the solution of problem ([I]) - ([5]) the following apriori estimate takes place 
(see, e.g., H) 

/ n \ 1/2 

||^lk{grad)<MK^||/,f • (6) 

\a=0 / 

Unstructured computational grids are often used to solve approximately 
boundary problems of mathematical physics. In some cases it is not suitable 
to use the coordinate-wise form of the considering equations and boundary 
conditions in one or another coordinate system. It seems more reasonable to 
employ the invarinat formilation of the problem. In this case the considered 
problem is represented via operators of vector and tensor analysis. 

Equation ([1]) is the basic one for problems of continuum mechanics. Vec- 
tors 

a={ao}, b = a = l,...,n, 
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are connected with the convective transport. Symmetric tensor of second 
rank 

D = {(io^}, a = l,...,n, /3 = l,...,n 

in ([1]) corresponds to the diffusion transport. Coefficient c(x) can be treated 
as a reaction. In such an interpretation equation ([!]) is nothing buth the 
equation of convection-diffusion-reaction 

— div(Dgrad u) + div(a'u) + bgrad u + c(x)'u(x) = /(x), x G ri. (7) 

Boundary value problem ([I]) - ([S]) in form ([7]) does not employ explicitly 
Cartesian coordinates 

3. Delaunay triangulation and Voronoi diagrams 

Assume that the computational domain is a convex polyhedron Q, n = 
2, 3 with boundary dQ. In domain n = ^lUdil we consider the grid To, which 
consists of nodes xf , i = 1,2, . . . , Md, and the angles of the polyhedron fl 
are nodes. Let a; be a set of inner nodes and du is a set of boundary nodes, 
i.e., uj = ZJ r\ Q, du = u D dfl. 

Each node xf, i = 1,2,..., M^, connect a certain part of the computa- 
tional domain, namely, the Voronoi polyhedron or its part belonging to Q. 
The Voronoi polyhedron (polygon) for a separate node is a set of points lying 
closer to this node than to all the other ones: 

Vi = {^\^eQ, |x-xf I < |x-xf I, j = 1,2,... ,Md}, i = l,2,...,MD, 

where | ■ | is the Euclidean distance. Each vertex x^, k = 1,2,..., My of 
the Voronoi polyhedron is associated with the tetrahedron constructed by 
the appropriate nodes contacting the Voronoi polyhedrons. We will assume 
that all vertices of the Voronoi polyhedrons lie either inside the computa- 
tional domain fl or on its boundary dQ. These tetrahedrons determine the 
Delaunay triangulation — a dual triangulation to the Voronoi diagram. The 
Z)-grid in the domain Q is determined by the set of nodes (vertices of tetra- 
hedrons of the Delaunay triangulation) xf, i = 1,2, . . . , Md, the V-grid is 
defined by the set of nodes (vertices of polyhedron of the Voronoi diagram) 
x^, k = l,2,...,Mv. 

We mark a separate tetrahedron Dj. of the Delaunay triangulation. This 
tetrahedron is identified by the number k of the Voronoi polyhedronm vertex. 
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k = 1,2,..., My. The tetrahedrons D^, k = 1,2, ... , My cover the entire 
computational domain, so 

_ Mv 

n= U Dk, Dk = DkUdDk, DfcnD„ = 0, k ^ m, k, m = 1,2, . . . , My. 

k=l 

For common planes of the tetrahedron we use the notations 

dDkm = dDkOdDm, ky^m, k, m = 1,2, . . . , My. 

The boundary of the computational domain dQ consists of the planes of 
Delaunay tetrahedrons. Let 

dDo = dn, dDko = dDkOdDo, k = 1,2, My. 

We associate the Voronoi polyhedron Vi, i = 1,2, . . . , Mjy with the node 
of the main grid i. Thus, we have 

_ Mo 

n= u v„ v, = v,udv,, Vinv,=^, i^j, t, j = i,2,...,Md 

i=l 

and 

dv^j = dv,ndVj, i^j, i, j = 1,2,... ,Md. 

For each tetrahedron D^, k = 1,2, . . . , My, we define the set of neighbors 
yV^{k), having common planes with D^, i.e.. 



W^{k) = {m I dDk n dDm 7^ 0, m = 0, 1, . . . , My}, k = l,2,...,M^ 



V- 



In this case, m = means that the tetrahedron contacts the bound- 
ary. We define also the neighbors for each Voronoi polyhedron Vi, i = 
1,2,..., Mn: 

W^W = {j I dVndV.^dS, j = l,2,...,Mn}, i = 1,2,..., Mn. 

We assume that the introduced Delaunay triangulation and the Voronoi 
diagram are regular 2l(|. For the notations 

= diam(Dfe) — diameter Dk, 

Qk = sup {diam(5') | S — sphere in Dk], k = 1,2, . . . , My 
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the regularity condition of the Delaunay triangulation is 

-^<(7>0, A; = 1,2,..., My. 

Likewise, for the Voronoi diagram we have 

hf — diam(V^) — diameter Vj, 
= sup {diam(S') | S — sphere in V^}, 
-^<a>0, i = 1,2,...,Md, 



h = max{/i^ 



and 



meas(Dfc) = / dx, k — 1,2, ... , My, 

meas(V^i) = / dx, i = 1, 2, . . . , M^. 
Jvi 

We will approximate the scalar functions of the continuous argument by 
the scalar grid functions that arc defined in the nodes of the D-grid or in the 
nodes of the y-grid. We denote by Hd the set of grid functions defined on 
the D-grid 

Ho^{y{^) I y(x) = ^(xf ) = , i = 1, 2, . . . , Mz, }. 
For the functions ?/(x) e Hd, vanishing on the boundary du, we define 

= { y(x) I y(x) e Hd, y(x) = Q,^eduj}. 

We consider the scalar product and the norm for the scalar grid functions 
from Hd by 

{y,v)D ^^yfv^ meas{Vi}, \\y\\D = (y,?/)^^ 
1=1 

This scalar product and the norm are grid analogs of the scalar product and 
the L2(0)-norm for the scalar functions of the continuous argument. 
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To determine the vector field in the control volume, it is natural to use 
the components of the sought function normal to the corresponding planes 
of the control volume. Choosing the initial and the final node, we connect 
with each tetrahedron edge D^, k = 1,2, . . . , My or polyhedron edge Vi, i = 
1,2,..., Md, a vector — the directed edge. For Delaunay triangulation the 
normals to the planes are the directed edges of Voronoi diagram and vice 
versa. For the approximation of the vector functions thereby we can use 
projections of the vectors on the directed edges. We will further use exactly 
this variant with the description of the vector field in the control volume by 
means of vector projections on the edges of the control volume. 

We will orient the Delaunay triangulation edges by the unit vector 

eg = ej^, i = 1,2,...,Md, J G >V^(^), 

directed from the node with a smaller number to the node of a larger number. 
The vector function y(x) on the Delaunay triangulation is defined by the 
components 

yj = yeg, i^l,2,...,Mn, i e >V^(i), 
that are given in the middle of the edges 

xS = ^(xf + xf). 

For the Delaunay triangulation used and the Voronoi diagram we define the 
length of the edges in the following way: 

/g = |xf-xf|, 1,2,..., Mn, je>V^(i). 

We denote hj Hd the set of grid vector functions determined by the 
components y^, i = 1,2, . . . , Mjy, j G W^(i) that are given in the middle 
of the edges. In a similar way we denote by Hy the set of grid vector 
functions defined by the components t/j^, k — 1,2,..., My, m e yV^{k). 
If the tangential components of the grid vector functions y e Hd vanish on 
the boundary, we define 

Ul^iylyeHn, y(x)eg = 0, x = xge9a;, 
1 = 1,2,. ..,Md, je>V^(z)}. 
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Consider the scalar product and norm 

Mo 

(y, v)d = XI 1^ '"S 1^^ ~ ^Sl meas(ai^ij), ||y||D = (y,y; 
in Hd- 

4. Gradient and divergence approximations 



,1/2 

'd ■ 



Let us present the main results (see [lOj for deatails) connected with 
approximation of gradient and divergence operators on the basis of the De- 
launay triangulation and Voronoi diagram. 

The set of grid functions Hd can be the domain of definition of the grid 
gradient operator. We denote the grid gradient operator by gradj). Taking 
into account the chosen edge orientation, at the points xf^ we set 



(grad^jy)g = 77(2,j) , i = 1,2, Mn, jeW'^it) 



i.e., the range of values of the operator grad^j : Hd — ?■ H^) is the set of vector 
grid functions Hd- In ([8]), we use the following notation: 



1, j > h 
-1, J < i. 



For the truncation error of the grid operator grad^ we have 

(grad^?/)(x) = (gradz/)(x) + g(x), g = 0{h^), x = xg, (9) 

1 = 1,2,. ..,Md, jeW^(2) 

provided that w(x) is a sufficiently smooth functions. 

For the Voronoi polyhedron this equality is written in the following form: 

divu dx = N / (u ■ n^) <^x, (10) 



where n]^ is the normal to the edge dVij outside with respect to V^. To 
construct the grid operator div/j : H^) — )■ Ho we use the elementary formulas 
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of integration for the left- and right-hand sides f llOl) . This leads to the grid 
analogs of the operator div in the form of 

(divDy)f = Kttt Y1 ("^S ■ e?)?/*?"^eas(ayij), i = l,2,...,MD. 

(11) 

We have 

(divDu)(x) = (divu)(x) + g{x.) + (div£,q)(x), 

g = 0{h), q = 0{h), x = xf, I = 1,2,..., Mn, (12) 

for the truncation error. The truncation error for the grid divergence oper- 
ators divD equals to C(l). However, there does exist a special divergence 
expression of the truncation error saving the situation in the case of approx- 
imation of problems of mathematical physics. 

Let us highlight the adjoint property of the considered grid operators of 
gradient and divergence. Taking into account the above notations we have 

{v, divD y)D + (grad^, v, y)/, = 
E E (4-eg)^^t±^^gnieas(ai-.,). (13) 
From (fT3|) it follows that 

div^ = -grad^ (14) 
on the set of grid functions v G i?^ and y E Hd- 



5. Approximation of the convection-diffusion-reaction problem 

In domain Q we consider the regular D-grid and y-grid. Taking into 
account boundary condition ([5]), we will find the approximate solution of 
problem ([5]), ([7]) as grid function y(x) G i/^. 

A particular discussion should be given to approximation of convective, 
and moreover, diffisive terms in equation ([7]). 

We consider approximations of the convective transport operator both in 
the divergent and in non-divergent (characteristic) form: 

Ci(v)'u = vgradw, C2(v)'u = div(vM) 

at specified vector field v(x), x G ^2. 
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The tangential components are approximated at the edge midpoints of 
the Delaunay triangulation as follows 



For scalar functions it is natural to introduce 

2/(xg) = ^(l/(xf) + 2/(xf)). 

The direct approximation of the convective transport operator in the diver- 
gent form results in 

Cf (v^)y = div^(v^?y) = ^— 1— Yl (ng-eg)^g(yf+yf)meas(9\/.,). 

(15) 

As for the convective transport operator in the non-divergent form, its ap- 
proximation is performed (l2| using the grid analog of the following equality 

Ci(v)m = C2(v)m — wdivv. 

At 

Cf(v> = Cf(v>-2/div^v^ 

we have 

(cf(v^)y)f = ^^^^;^ ^^^K-^S)-Siyf -y^)^-^<9v.,). (le) 

For the trancation error similarly to f|T2l) we have 

C^{v^)u = Ca(y)u + 5f(x) + (divDq)(x), a = 1, 2, 

g = 0{h), q = 0{h), x = xf, i = 1,2,...,Md. (17) 
Let us construct a grid analog of the diffusion transport operator 

Vu = div(D grad u), x G 17. 

We put formally 

P^?/ = divz5(D^grad^y). (18) 
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Figure 1: Two-dimensional case 



Above (see ([H])) we presented the grid approximation for the projection of 
the gradient operator onto edges of the Delaunay triangulation. It is enough 
for approximation of the diffusion operator in isotropic media (with a scalar 
diffusion coefficient). For general elliptic equations with mixed derivatives ([1]) 
(with a tensor diffusion coefficient) it is necessary to approximate the gradient 
not only along the edges of triangulation. To investigate this problem, let us 
start from the 2D case (n = 2, FigH]). 

Here i,j,k,m — nodes of the Dalaunay triangulation, i,j,k i,j,m — 
two adjacent triangles of this triangulation with common edge p q — 

nodes of the Voronoi diagram. For nodes of the Voronoi diagram there are 
used notations = x^^, and x^ = x^^. Highlighted in FigJT] edges of the 
Delaunay triangulation and Voronoi diagram are orthogonal to one another. 
Introduce at the point of their intersection (at point x = x^) local coordinate 
system (si,S2) with direction unit vectors 61,62). 

Divergence of a vector field at the discrete level (see (fTTj) ) is defined 
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Figure 2: Tree-dimensional case 



using the tangential component. To approximate the diffusion transport, it 
is necessary to evaluate the diffusion tensor ( ITSll . Assume for the introduced 
local coordinate system 

Coefficients {d^p)ij, a, (3 = 1,2 are calculated via specified smooth enough 
coefficients dapi^x), Oi,j5 = 1,2 of equation ([1]). If local coordinate system 
(si,S2) coincides with (xi,X2) we have, for instance, 

ida^)ij = dap{x?j), a, /3 = 1, 2. 

In the general case the recalculation rule for a tensor of second order should 
be taken into account when going to a new orthogonal coordinate system. 

Using the above notations in local coordinate system (si, S2) the tangen- 
tial component of the flow can be written as 

(Dgrad^y)g = (rffjij(grad^ + ((if2)ii(grad^ 2/)2. (19) 
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The tangential component of the gradient seems hke this 



(grad^ y)i = (grad^y)g = ^1—^. (20) 



For the normal component we have 

(grad^?/)2 = ^V^' (21) 

where = |x^^ — x^^| — the edge length of the Voronoi polygon (polyhe- 
dron) . 

Evaluation of the normal component of the gradient via (12T|) is based on 
interpolaton of values of a scalar function defined at the Delaunay triangula- 
tion nodes to the Voronoi polygon nodes. Let Ipij = |Xp — x^| be the distance 
from node p of the Voronoi partitioning up to edge {i, j) of the Delaunay tri- 
angulation so that = dptj + dgij. The linear interpolation with accuracy 
0{h?) over triangle Dp results in 



y^p = rr..l(Ti\ iy^^^dp^^ + y?^%dp^>^ + yfi^dp^k) ■ (22) 

meas^ i-jp j 



Taking into account ( l20|) - ( !22l) . the next representation for the truncation 
error of the flow follows from (fT9l) 



(Dgrad^«)g = (Dgrad w)(xg) + 0{h). (23) 
Approximation ([T]) leads to the grid equation 

- divz5(D^ grad^ y) + Cf (a^)y + Cf (b^)|/ + c^y = 

x = xf, z = 1,2,...,Mb, (24) 

where, for instance, c-^(xP) = c'^(xf). Grid equation ( 124|) approximates 
equation ([1]) with error 

^(x) = (7(x) + (div,,q)(x), g = 0{h), q = 0{h). (25) 

In the 3D case approximations of the diffusion transport operator are 
constructed in a similar way. Let us consider (see Figj2]) two adjacent tetra- 
hedrons of the Delaunay triangulation — i,j,k,m and i,j,k,n. For the 
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corresponding nodes of the Voronoi diagram (vortices of the Voronoi poly- 
gon) we again use notations and x^. The hne connecting these vortices 
is perpendicular to triangle k and intersects it at point r (Figj2]). Point s 
is the midpoint of edge Stright line (r, s) is perpendicular to both edge 

of the Delaunay triangulation and edge {p,q) of the Voronoi diagram. 
On the basis of these three lines {p,q) and (r, s) highlighted in Figl2] 

we construct local coordinate system (51,52,53) with direction unit vectors 
61,62,63). 

Similarly to the 2D case we can evaluate now the diffusion tensor in the 
local coordinate system 

({dn)ij {dr2)ij (^13)^ \ 
(dpj (dgh (dgh ■ 
('^3l)ij {d^2)ij {d^^)ij J 

For the smooth enough coefficients c?q,^(x), a,j3 = 1, 2, 3 of equation ([T]) and 
in case of coincidence of local coordinate system (51,52,53) with {xi,X2,X3) 
it is possible to take 

idal3)ij = daisi^fj), a, ^=1,2, 3. 

In local coordinate system (51, 52, 53) the flow component along edge (z, j) 
of the Delaunay triangulation is equal to 

(D grad^ y)g = 

{dii)ij{gr&djjy)i + (c?f2)ii(grad^ y)2 + {d^^)ij{gTad^y)3. (26) 

The tangential component of gradient (grad^y)i is calculated in accordance 
with (I20!) whereas for component (grad^, y)2 (along the edge of the Voronoi 
diagram) approximation (12T|) is used. For the third term in (126|) we employ 

(grad^?/)3 = ^V^, ^27) 

where Irs — distance between points r and 5. Values of scalar function yr, tjs 
in (127|) are approximated with the second order via ( l22|) -type formulas using 
values defined at nodes of the Delaunay triangulation. 

So, representation f l23p for the truncation error of the flow is valid in the 
3D case, too. Hence we aso obtain grid equation flMj) with representation 
f l25|) for the truncation error. On the basis of this divergent representation 
for the error we can prove convergence of the discrete solution to the exact 
one [lO| in the grid analog of space 'H(grad) with the first order with respect 
to h. 
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